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The exact nonequilibrium steady state solution of the nonlinear Boltzmann equation for a driven 
inelastic Maxwell model was obtained by Ben-Nairn and Krapivsky [Phys. Rev. E 61, R5 (2000)] 
in the form of an infinite product for the Fourier transform of the distribution function /(c). In 
this paper we have inverted the Fourier transform to express /(c) in the form of an infinite series of 
exponentially decaying terms. The dominant high energy tail is exponential, /(c) ~ Aq exp(— a\c\), 
where a = 2/yl — a 2 and the amplitude Ao is given in terms of a converging sum. This is explicitly 
shown in the totally inelastic limit (a — > 0) and in the quasi-elastic limit (a — > 1). In the latter case, 
the distribution is dominated by a Maxwellian for a very wide range of velocities, but a crossover 
from a Maxwellian to an exponential high energy tail exists for velocities |c — co| ~ around a 

crossover velocity Co — \nq~~ 1 / yfq, where q = (1 — a)/2 <C 1. In this crossover region the distribution 
function is extremely small, ln/(co) — q~ Inq. 

PACS numbers: 45.70.-n, 05.20.Dd, 51.10.+y 



I. INTRODUCTION 

In kinetic theory there is a long standing interest in 
overpopulated high energy tails of velocity distribution 
functions Q because of chemical reactions and other ac- 
tivated processes that occur only at energies far above 
thermal. This interest has been considerably increased 
in the past 10 years because of research in granular fluids 
with dissipative or inelastic interactions. The velocity 
distributions in fluidized systems have been studied the- 
oreticaUy 110, S S S 0, 0, and measured in Monte 
Carlo |3.l9lll0| and molecular dynamics simulations fTTj| . 
and in numerous laboratory experiments |12| |. 

Very recently a revival in this field occurred when 
Baldassarri et al. 0, 0] discovered an exact scal- 
ing solution - with an algebraic high energy tail - of 
the non-linear Boltzmann equation for an inelastic one- 
dimensional freely cooling gas (without energy input) 
with a collision frequency independent of the energy 
of the colliding particles. This model, called Inelastic 
Maxwell Model (IMM), was introduced by Ben-Nairn 
and Krapivsky 15]. It is in fact an inelastic modifica- 
tion of Ulam's stochastic model to illustrate the velocity 
relaxation of elastic one-dimensional point particles to- 
wards a Maxwellian [ljj} . A three-dimensional version of 
it has been constructed by Bobylev et al. 0, 0] . For 
a recent review on inelastic Maxwell models, see Refs. 

EH. 

Baldassarri et al. have demonstrated the importance of 
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this type of solutions in 0| with the help of Monte Carlo 
simulations of the nonlinear Boltzmann equation for 1- 
dimensional and 2-dimensional IMM's. It appeared that 
the solution, F(y,t), for large classes of initial distribu- 
tions F(v, 0) (e.g. uniform or Gaussian) and for all values 
of the inelasticity could be collapsed for large times on 
a scaling form VQ d (t)f(v/vo(t)), where Vo(t) = (v 2 ) 1 / 2 
is the r.m.s. velocity. In one dimension the scaling 
form was given by /(c) = (2/7r)(l + c 2 )~ 2 , which has 
a heavily overpopulated algebraic tail ~ c~ 4 when com- 
pared to a Maxwellian. In two dimensions the solu- 
tions also approached a scaling form with an algebraic 
tail, /(c) ~ c~ d ~ a with an exponent a(q) that depends 
on the degree of inelasticity q — ~(1 — a), where a 
is the coefficient of restitution. Soon after, Ben-Nairn 
and Krapivsky pl[. and Ernst and Brito (2^ obtained 
asymptotic solutions with algebraic tails for the velocity 
distribution in d-dimensional freely cooling IMM's from 
self-consistently determined solutions of the Boltzmann 
equation. Using methods previously developed for the 
inelastic hard sphere case, the asymptotic solutions were 
also extended to non-equilibrium steady states (NESS) 
in d-dimensional systems driven by Gaussian white noise 
and other thermostats [H |2^] . There the tails exhibited 
over-populations of exponential type, ~ exp(— a\c\), for 
all (i-dimensional IMM's [23j. For inelastic hard spheres, 
which is the prototypical model for granular gases, the 
velocity distribution function shows an overpopulated ex- 
ponential tail in free cooling 0, 0. @, and a stretched 
exponential tail ~ exp(— et|c| 3 ' 2 ), when driven by white 
noise H0G3- 

For the case of d-dimensional free IMM's the approach 
of F(v, t) to a scaling form with an algebraic tail has also 
been rigorously proven, for initial distributions in the C\ 
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function space, satisfying the physical requirements of 
finite mass and energy, i.e. J dv{l, v 2 }F(v, 0) < oo |24| . 

What about exact and/or more explicit results for the 
distribution function in the onc-dimcnsional IMM, driven 
by Gaussian white noise? The exact solution of the non- 
linear Boltzmann equation for this case is given in the 
form of an infinite product for the Fourier transform of 
the distribution function |l5j . Nienhuis and van der Hart 
[2q made an extensive numerical analysis of this solution, 
and demonstrated exponential decay, in agreement with 
the predictions of Ref. 23;]. More numerical evidence 
for exponential high energy tails in the one-dimensional 
driven IMM was given recently by Marconi and Puglisi 
pfij , and by Antal et al. • I n a recent paper [23] , Ben- 
Nairn and Krapivsky have also used the Fourier trans- 
form method to show that the high energy tail is expo- 
nential for any inelasticity, but with an amplitude that 
diverges in the quasi-elastic limit. On the other hand, 
the problem of determining for what range of velocities 
the exponential tail actually applies remains open. This 
is one of the points addressed in this paper. 

The plan of the paper is as follows. In the remain- 
der of this Section we present the nonlinear Boltzmann 
equation for the velocity distribution function F(v) or 
/(c), driven by Gaussian white noise, and we discuss 
qualitatively the physical properties of the model in dif- 
ferent limiting cases. In Section [H] the exact solution 
<f>(k) = J dce~ f(c) of the Fourier transformed Boltz- 
mann equation in the NESS is presented in the form of an 
infinite product, and its large- and small- k properties are 
analyzed. In Section lTTTl we determine the inverse Fourier 
transform, /(c), in the form of an infinite series of expo- 
nentially decaying terms. In the limit of totally inelastic 
collisions (a — * 0) substantial simplifications occur. The 
rather singular quasi-elastic limit (a — > 1) is studied in 
Section llVl where also the crossover from Maxwellian to 
exponential decay is analyzed. We end with some com- 
ments in Section[V] and some technical details are moved 
to Appendices lAl and iBl 

Before concluding this introduction we present the 
Boltzmann equation for the one-dimensional inelastic 
Maxwell model (IMM) |l5j. driven by Gaussian white 
noise, and we discuss some of its important properties. 
The time evolution of a spatially homogeneous isotropic 
velocity distribution function F(v,t) — F(\v\,t) is de- 
scribed by the nonlinear Boltzmann equation , 

- D*$* = Jd Vl [±F(v»)F(v>{) - F(v)F(vi)] 
= -F(v) + \JduF{u)F = I(v\F). (1.1) 

All velocity integrations extend over the interval 
(— oo, +oo). The diffusion term represents the (heating) 
effect of the Gaussian white noise with noise strength 
D. The nonlinear collision term represents the inelas- 
tic collisions, where v" = v — |(1 + a^ r )(v — Vi) and 
v'l = vi + |(1 + a _1 )(-y — vi) denote restituting veloc- 
ities. Here a = 2p — 1 = 1 — 2q with < a < 1 is 
the coefficient of restitution. The mass is normalized as 



J dvF(v) = 1, and the mean square velocity or tempera- 
ture as (v 2 )(t) = J dvv 2 F(v) = v 2 (t) . The rate equation, 

d t (v 2 ) = 2D-2pq(v 2 ), (1.2) 

obtained from Ijl.lf) . describes the approach to the non- 
equilibrium steady state (NESS) with width (v 2 ) = 
D/pq, where the heating rate D caused by the random 
forces is balanced by the loss rate, pq(v 2 ) = j(l—a 2 )(v 2 ) , 
caused by the inelastic collisions. 

To understand the physical processes involved, we first 
discuss in a qualitative way the relevant limiting cases. 
Without the heating term (D = 0), Eq. (|1.1|) reduces to 
the freely cooling IMM, whose exact solution has been 
discussed in Refs. hjhj. If one takes in addition the 
elastic limit (a — > 1 or q — > 0), the collision laws reduce 
in the one- dimensional case to v" = Vi,v" = v, i.e. an 
exchange of particle labels, the collision term vanishes 
identically, every F(v, t) — F(v) is a solution, there is no 
randomization or relaxation of the velocity distribution 
through collisions, and the model becomes trivial at the 
Boltzmann level of description, whereas the distribution 
function in the presence of infinitesimal dissipation (a — > 
1) approaches a Maxwellian. 

If we turn on the noise (D ^ 0) at vanishing dissipation 
(q = 0), the exact solution of (| 1 . 1 1) in Fourier represen- 
tation is F(k,t) = ex.p(—Dk 2 t) F(k,0), and the granular 
temperature, v 2 (t) = v 2 (0) + 2Dt, increases linearly with 
time. With stochastic heating and dissipation (even in 
infinitesimal amounts) the system reaches a NESS, and 
it is the goal of this paper to determine the NESS distri- 
bution function. 

To expose the universality of this NESS it is convenient 
to measure the velocities, c = v/vq(oo), in units of its 
typical size v (oo), i.e. the r.m.s. velocity or width of the 
velocity distribution uo(oo), 

F(v,oo)=vv 1 {oc)f(v/v {oo)), (1.3) 

which obeys the normalizations J dc {1, c 2 }/(c) = jl, 1}. 
Different normalizations have been used as well |28| . 

The rescaled velocity distribution in the NESS is then 
the solution of the scaling equation, 

J( C |/) = --A-/"( C ) = -pqf'(c), (1.4) 
v 0\°°) 

where primes denote c-derivatives. The first equality may 
suggest that /(c) may depend on the noise strength D 
and possibly on the initial distribution via i>o(oo). By 
eliminating wo(oo) with the help of 11. 2|) in the NESS we 
have shown that the scaling form of the distribution func- 
tion /(c) is a universal function, that does not depend on 
the strength D of this thermostat, nor on any property 
of the initial distribution. It only depends on the type of 
thermostat used. 
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II. FOURIER TRANSFORM OF IMM 
BOLTZMANN EQUATION 

The nonlinear Boltzmann equation for characteristic 
function, 4>(k) = J dc e~ lkc f (c) , is obtained by Fourier 
transformation of (|1.4f> with the result, 

(1 + pqk 2 )(j}{k) = (j){pk)(j){qk). 



In particular, C2 = (c 2 ) = 1. Since 1 — p 2n — q" n > 
0, it follows that all cumulants are positive, indicating 
already an overpopulation of the high energy tails. So far 
a summary of the results obtained in Ref. [l5[ . We note 
that the Stirling approximation shows that the cumulants 
at fixed a or q and n > e/(2y/pq) are rapidly diverging 

with increasing n, as C2n ~ 2y/ / n(2n^/pq / e) 2n . 

The exact solution 4>(k) in 1|2.3|1 has an infinite se- 
quence of poles of multiplicity v m g in the complex k- 
plane, all of which contribute to the amplitude of the 
asymptotic high energy tail of /(c). This makes a nu- 
merical inversion of <p(k) to obtain f(c) a bit tricky. To 
determine /(c) several authors pH |27| have performed 
numerical inversions of <fi(k), starting from the infinite 
product (|2.3[l or from the more convenient series expan- 
sion 12.4(1 . However the latter one is only convergent for 



■2,1 



(2.1) 



The simple structure of the equation for the Fourier 
transform (p(k) follows because the nonlinear collision op- 
erator for (in) elastic Maxwell models is a convolution in 
the velocity variables pj. Equation (|2.1|l is a nonlinear 
finite difference equation, that can be solved by itera- 
tion. A simple way to construct the exact solution is to 
introduce ip(k) = ln</>(fc), which satisfies 



rjj(k) = ip(pk) + ip(qk) — In (l + pqk 2 



The normalization of mass and energy implies that 
<f>(k) « 1 - \k 2 and ip{k) m ~\k 2 at small k. The 
solution to (|2.2(1 can be found iteratively starting from 
ipo(k) = — ln(l +pqk 2 ) and inserting ip n (k) on the right- 
hand-side of 11'. 21 to get ip n +i(k) on the left-hand-side. 



By taking the limit ip{k) 
obtains, 



limn^oo ip n (k), one finally 



(2-2) pqk 2 < 1. To facilitate such numerical procedures, we 



have derived an expansion in powers of k~ , convergent 
in the complementary region, pqk 2 > 1, of the complex 
/c-plane. This rather technical part is deferred to Ap- 
pendix^ The results can be found in l|A2(l . I|A4(1 . and 
ifATjl . 



III. HIGH ENERGY TAIL 



<p(k) 



00 m 



(2.3) 



m=0 e=o 



where v m t = (™). These solutions satisfy the required 
boundary conditions at k = 0. We further note that 
ip(k) = ip(k) — \\k\ with A an arbitrary complex number 
is also a solution of 12. 2|) . but in general does not sat- 
isfy the boundary conditions at small k. This property 
is a reflection of the Galilean invariance of the original 
Boltzmann equation. 

Equations (|2.3|) provide an exact representation in 
Fourier space of the solution of the Boltzmann equation 
The series (|2.3|) converges rapidly, even for large 
k. By expanding the logarithm in powers of k 2 and 
summing a geometric series we obtain, 



y, (-1)" (k 2 pqr 

n=l 



II 



1 — p 2n — q 



In ' 



(2.4) 



It converges for k 2 < 1/pq, and ip(k) has a branch point 
singularity at k 2 = —1/pq, as is apparent from 12.211 . 
Equation (|2.4|) allows one to get the cumulants Cm de- 
fined by 



4>(k) 



OO 



with the result 



c 2 



n—l v ' 



(2n)! (pq) n 



1 — p 2n — q 2 



(2.5) 



(2.6) 



On account of i|2.3|) . the characteristic function cj>(k) 
can be written as 



2 \-"m< 

me) 



(3.1) 



m=oe=o 



where k m i = ap~ e q~( m ~ e ) with a = 1/ \fpq~. Thus 4>{k) 
has poles at k — ±ik m i with multiplicity v m i. The ve- 
locity distribution, 



(3.2) 



can then be obtained by contour integration. As /(c) is 
an even function, we only need to evaluate the integral in 
(|3.2|) for c > 0. The replacement c — !■ |c| gives then the 
result for all c. By closing the contour through an infinite 
upper half-circle and applying the residue theorem we 
obtain 



m -EE 



00 m — 1 

k m e\c\ 



E \ c \ nA 



III < II • 



(3.3) 



m=0 1=0 



n=0 



where 
with 



1 ^mi 



lim 



d_ 

n\(y m g — 1 — n)\ k^ik mt \dk 
x(k + ik mi )-^ e me (k), 



(3.4) 



w=nn (i+k 2 /k 2 m , t ) 



m'=0 l'=§ 



(3.5) 
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Note that the factor labeled (m',£') — (m,£) is absent. 
The dominant terms in 13.311 for large |c| correspond to 
the smallest values of k m i- The two smallest ones are 
fcoo = a and kn — a /p. Consequently, the leading and 
sub-leading terms are 



f(c)^A a e- a ^+A ie 



,-a\c\/p 



where 



A,, 



Anno = {a/2p n )(/)nn(ia/p n ). 



(3.6) 



(3.7) 



We calculate the first two explicitly, i.e. 



An = 



A, 



2 6XP 



E P + 1 
m (1 - p 2m - q 2m 

ra=\ v 1 H 



(3.8) 



-ap 



2(1 -p 2 )(p-q) 



x exp 



- p -2m(p2m + g 2m)2 

m (1 



.m— 1 



P 



2m n 2m 



q 2m ) 



(3.9) 



In the last equalities we have followed steps similar to 
those used to obtain (|2*1J) from i(2~3|) . The results (|3~3|> - 
(|3.9|) exhibit the full analytic structure of the dominant 
and subdominant high energy tails of the velocity dis- 
tribution in the NESS, as already demonstrated numer- 
ically for the one-dimensional case in Refs. I2H |27|. 
and derived in |2jj for d-dimensional IMM's on the ba- 
sis of self-consistent solutions. Moreover we have ob- 
tained here explicit expressions for the amplitudes Aq 
and Ai in the form of sums that are rapidly converging 
when q is not too small. The coefficients Aq = Aqoq and 
Ai = Auo are shown in Fig. ^ as functions of a, where 
Ano oc l/(p — q) — 1/a diverges according to l|3.9|l . The 
next term to those explicitly given in 1(3. 6|) corresponds 
either to fc 2 2 = a/p 2 if p 2 > q (i.e., if a > V5-2 ~ 0.236) 
or to fcio = a/q if p 2 < q. Note that the amplitude Aioo 
of exp(— fcio|c|) can be obtained from Ai in (|3.9|l by in- 
terchanging p <-» q. Figure [21 compares the asymptotic 
form /(c) ps yloe _a l c l with the function /(c) obtained by 
numerically inverting <p(k) for a — and a = 0.5. We ob- 
serve that the asymptotic behavior is reached for a\c\ > 4 
if a = and for a\c\ > 8 if a = 0.5. As a = l/s/pq this 
corresponds to velocities far above the r.m.s. velocity. 

There are two interesting limiting cases: the quasi- 
elastic limit (a — ► l,q — ► 0), and the totally inelastic 

limit (a — > 0,p — * ^ + ,q — > | ). We start with the 
latter. In the totally inelastic limit (a — > 0) the subdom- 
inant terms Aiioe _a ' c '/ p and Aiooe _a ' c '/ 9 become equally 
important, i.e. the single poles in 1)3.1(1 at kn = a/p 
and fcio = a/q coalesce, and (13.611 no longer describes 
the subdominant large-c behavior correctly. Moreover 
Auo — — Aioo oc 1/a, as can be seen in Fig.^for ^uo- In 
fact, the poles k m i — > k m = 2 m a coalesce for all £, some 
of the coefficients A min diverge, e.g. A mm0 cx (1/a) 2 _1 , 
and the expansion makes no sense anymore. So, we an- 




0.0 0.2 0.4 0.6 0.8 1.0 



a 



FIG. 1: Logarithmic plot of the amplitudes Aq = Aooo (solid 
line) and —A\ = —Auo (dashed line) as functions of the 
coefficient of restitution. The arrow indicates the value An ~ 
2.958389 at a = 0. The dotted lines represent the asymptotic 
form 14. 1H for small q. 




a Id 



FIG. 2: Logarithmic plot of /(c) versus a\c\ for a = and 
a = 0.5. The dotted lines are the asymptotic forms /(c) « 
Aoe _a ' c ' at a — and a = 0.5, with Aq obtained from l|M.8|l . 



alyze the case a — separately. In this case the charac- 
teristic function is according to l|2.3|) . 



m= n (i+fcv^r™, (3.io) 

m=0 

where = 2 m and fc m = 2 m a with a = l/^/pq = 2. 
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Then the distribution function is 

oo u 77l — 1 

/(c) = £ e- fc ™l c l £ MM™, 

m— n— 

where the residues or amplitudes are given by, 



(3.11) 



A 



i n+1 k 2 ^ / 5 

hm 



n\(y m — 1 — n)! k^ik m \dk 
x(k + ik m )- Um m (k), 



and 4> n (k) is defined as, 



t (fc)= rj (i + fc 2 /fc. 



2 /J.2 \-' / m(l-5nro) 



(3.12) 



(3.13) 



m=0 



For large |c| the distribution function becomes, 

/(c) « A 00 e- 2 I C I + (A 10 + A u |c|)e- 4 I C I + ■ • • . (3.14) 

To calculate the amplitudes of the dominant terms we 
derive from l|3.13|l . 

oo 

ln^ (fc) = - 2 m ln(l + 2- 2m fc 2 /a 2 ) 

771—1 

(-1)" (fc/a) 2n 



E 



n 2 2 ™" 1 - 1 



ln0i(fc) = -ln(l + /c 2 /a 2 ) 



+2E 



(-1)™ (fc/2a) 2 " 



n=l 



2 2n-l _ 1 • 



(3.15) 



The definitions <|3~T21l - (|3~T5|> with fc m = 2 m a(a = 2) yield 
then, 

Aqo = \a4> (ia) = e s ° ~ 2.958389 

An = a 2 0i(2ia) = -|A§ ~ -11.669422 

Aiq — ^a<j>i(2ia) — ia<f)\{2ia) 

= | (Si -§) Al Q ~ 3.138267, (3.16) 

where we have used the rapidly converging sums, 

5 = V-(2 2n " 1 -l)- 1 - 1.084645 

n—l 

oo 

51 = J](2 2n - 1 - I)" 1 ~ 1.185597. (3.17) 



n=l 



In fact, the results (|3.14l) could have been derived directly 
from l|3.3|) - (|3.9|) after lengthy calculations, by expanding 
j4ho and Aiqq in powers of a, with the result, 

A ls0 = (-l) s+ Mn/(8a) + \A W + 0(a) (3.18) 

with s = 0,1. Insertion of these results in l|3.3|l yields 
(|3.14|) . The limit a — > 1 is discussed in the next Section. 



IV. QUASI-ELASTIC LIMIT 

As already mentioned in the introduction, the velocity 
distribution approaches a NESS, even in the presence of 
an infinitesimal dissipation (a — > l,q — > 0), balanced 
by a ditto amount of stochastic heating. This limit is 
referred to as the quasi-elastic limit. For the rescaled 
functions /(c) and <fi(k) it simply refers to the limit q — > 0. 

Once we have first taken the large |c| limit at fixed 
a < 1 — as has been done in the previous Section — we 
can next take the quasi-elastic limit a — > 1. When the 
limits are taken in that order, the asymptotic behavior 
is still of the form e~ a ' c ', where the decay constants are 
kmm = a/p m — > a, and the amplitudes may diverge. On 
the other hand, if the limits are taken in the reverse order, 
first a — > 1 at fixed |c| and next |c| — > oo, the behavior 
is in general totally different. 

First consider the second case, and observe that 
ijj(k) in H2.4fl has at small q the form ip(k) = 
— \k 2 + Y^=2 a ^n{q)k 2n with rapidly decreasing coeffi- 
cients a 2n ~ (-l) n g n_1 (l - ig)/(2n 2 ) for n > 2. Conse- 
quently <^>(fc) = e^C ) can be expanded as, 



1 



^2^2n(q)k 2 



n=2 



(4.1) 



where the relation between et2 n and fx^n is the same as 
between cumulants and moments after setting a,2 = \i2 — 
0. The coefficients [i2n are to dominant order in q 2 given 
by 



,i 4 = a 4 = \q{\ -\q- \q 2 ) + 0(q 4 ) 
= a 6 = -^g 2 (l-^)+0(9^ 

/'8 



\a 2 



£il2 



18^ 

4^ a s = i|gg 2 (l + 3 g ) + 0(g 4 ) 

a 4 a 6 = - is? 3 + 0(q 4 ) 
l 
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U2 = s fe y « 3 + o(? 4 ) 



(4.2) 



and in general /U4n-2 ~ M4n ~ C(?") for n > 2. The 
series above can be Fourier inverted term-wise, using the 
following relation, 



J — OO 27T 



ifc 2 fc 2„ 



= (-l)™(A) 2n ex p ( -I c 2)/^ 
He 2n (c)f (c) = 2"n!it 1/2) (ic 2 )/ (c),(4.3) 



= (-1) 



where /o(c) = exp(— ic 2 )/\/27r is the Maxwellian. In 
the last two equalities Rodrigues' formula for the Her- 
mite polynomials has been used, as well as their relation 
to the generalized Laguerre or Sonine polynomials (see 
Ref. |2<j, Eqs. (22.11.88), (22.5.18) and (22.5.40)). The 
resulting Sonine polynomial expansion of the velocity dis- 
tribution in the NESS reads then, 



/(c) - hie) 



l + J2(-l) n V2n(q)He 2n (c) 



(4.4) 



Similar expansions of the NESS-distribution function in 
low order Hermite or Sonine polynomials have also been 
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derived for inelastic hard spheres in d-dimensions [j| , and 
for a 3-dimensional IMM in |lSj . 

Next we consider the case, where first |c| — ► oo at finite 
a < 1, and ne:r£ a — ► 1, or <7 — * 0. The large-c behavior at 
fixed a has already been discussed in (|3.3I) - (|3.9I) . and we 
observe that the terms in Eq. (|3.3(l at large c, associated 
with all poles of the form k n e = a/p e q n ~ e (£ < n) decay 
rapidly as q — > 0, and only poles with k nn = a/p n need 
to be considered: 



f(c) = ]T A n e- k ™ 



(4.5) 



n=0 



We will analyze the behavior of the associated amplitudes 
A n by combining l|3.7|) with l|3.5|l . i.e. 

In 4, = ln(a/2p n ) 

- E7n=o Etc (7) (! - in [i - v<">-«] 

= + Si 2) + Bi 3) + ln(a/2p") , (4.6) 



where 



_ -2(n-<)^2(m- 



Tt — 1 III / \ 

m=0 f=0 V / 
n— 1 / \ 

B i 2) - -E(>[i-(^) 2M) 

B (3) = _ 



£=0 v 7 
oo m 



m=n+l 1=0 



(4.7) 



Now we take the limit q — > at finite n and retain terms 
to order q. The dominant small-g contribution to Bn^ 
comes from £ = m, i.e. 



^ln(l-p- 



2m \ 



o(q) 



In [(-2g) n n!] - k(n + 2)q + o(q), (4. 



where we used the relation 1 — l/p 2m ~ —2mq[l + (m + 
\)q\, and o(q k ) denotes terms which are negligible with 

respect to q k . Note that the complex number B„ is 
only determined modulo {2iri}, but exp is single- 

valued. Furthermore, we observe that = 0(nq 2 ). 

The analysis of B^ in igTJ) is more involved and given 
in Appendix FBI The result is, 

2 

4 3) - ^- + 5 ln «-^o+i (n+ § - g+o( ? ), (4.9) 



where 



3 7T 2 

*T = - + — - \ In 2- i?~ 0.733598. 



(4.10) 



Combining the small-g results l|4.8(l and i|4.9[l for B^ 
and i?i 3) with (|4.fijl yields for A n , 



(-1)" 

2n!(2<;)>« eX P 



exp + 0(g)] 

- K - \n{n - l)q + K x q + o(q) 



>(3) 



where 



25 
24 



144 



0.9731278. 



(4.11) 



(4.12) 



To describe the crossover between the two different lim- 
iting behaviors, i.e. I|4.4|) with first q — > 0, next c — > oo, 
and (|4.5(1 . (14. lip with first c — > oo, next g we need 
to couple these limits, which will be discussed next. 

By an extension of the steps followed in Appendix 
IBI it can be verified that the terms denoted by o(q) in 
Eq. (|4.11|) have the form n kl q k2 with ki < k 2 + 1 and 
k 2 > 2. Therefore, those terms can be neglected against 
the terms of order q if n <C q~ x . 

The ratio R(c) between the distribution function /(c) 
in (|4.5|l and its asymptotic high energy form Aoe~ a ' c ', 
defines a crossover function 



R(c) = /(c) M e- a l c l = j^6„r n , 



(4.13) 



n=0 



where r„ and 6 n = A n /Ao follow from l|4.5(l and (|4.11|l 



exp [-a\c\(p- n - 1)] 



(-1)" 
n,!(2q) n 



exp [-|n(n - 1)? + o(n 2 g)] . (4.14) 



Here wc have written o(q) — > o(n 2 q) to emphasize the 
fact that Eq. (|4.14|l remains valid if n <C q^ 1 ■ So there 
is a crossover behavior in R(c) from a large-c behavior 
of 0(exp[— c 2 /2]) ~ in the small-q Sonine polynomial 
expansion l|3.4|l . to the the small-g behavior of R(c) of 
0(1) in H4.13fl . The transition region is characterized 
by a crossover velocity Co such that R(cq) ~ h. The 
interesting questions are, how does Co scale with q in the 
quasi-elastic limit, and what is the width of the crossover 
region? To address these questions, note that the series 
(|4.13(l converges for all velocities and the signs of the 
terms are alternating. Therefore, when breaking off the 
infinite sum at n = TV, the maximum error is |6/v+i |rjv+i: 

AT 

R{c) = ^„ r „+AW( c ) 

= i?W(c) + AW(c), |AW(c)| < \b N+1 \r N+1 . 

(4.15) 

This suggests that the pure exponential high energy tail 
Aoe~ a l c l qualitatively describe the large-c behavior of 
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/(c) if 



l&ilri = 



-a\c\q/p 

2q 



~2q~ 



< 



(4.16) 



Of course the bound ^ may be replaced by any num- 
ber of the order of 1 in this estimate. Equation (|4.16() 
implies that w = \c\ Jqj \rtq~ 1 > 1. Therefore, we can 
estimate the crossover velocity to be Co = (In q -1 )/ yfq or, 
equivalently, Wq = 1. To confirm this and get a closed 
form for the crossover function R(c), consider a value 
of w in the range 0.5 < w < 1 and take N = j3q w ~ 1 , 
where > 1. In that case, TV > 1 but N 2 q « 1, so 
that r„ ~ q nw and £>„ ~ {-l) n / n\{2q) n for n < N, and 
|6jv+ikiV+i - (2/3/e)- Ar /2/3V27r^. Therefore, with this 
choice of N, 



N 



R (N) (c) - £ 

n=Q 

|AW(c)| < 



n / „w — l 



2^ 

2/3V2^/V V e 



-JV 



(4.17) 



If j3 > e/2 ~ 1.36 then AW(c) < 1 and R(c) can be 
approximated by R^ N '(c). By the same arguments, the 
upper limit in the summation of l|4.17|l can be replaced 
by infinity. The choice N = Pq w ~ 1 is justified by the fact 
that for w < 1 the term |6 n |r„ reaches a high maximum 
value \b no \r no ~ exp(n + |)/V27m at n ~ ^(g"' -1 - 
1) and then decays rapidly. If w > 1, however, |6„|r„ 
decreases monotonically and thus A^ N '(c) <C 1 for any 
choice of N. In conclusion, the crossover function for 
w > 0.5 in the quasi-elastic limit becomes 



R(c) ~ exp (-<f 



72) 



Icl^/lng- 1 . (4.18) 



At w = 1 we have i?(c = Co) — 1/ -y/e — 0.6, thus confirm- 
ing the estimate of the crossover velocity Co made below 
Eq. (|4.16() . Figure |3] represents the crossover function 
R(c) versus the scaled velocity w for q = 0.01, 0.001, 
and 0.0001. To measure the width of the crossover re- 
gion, let Wi and W2 denote the values of w at which 
R = 0.1 and R = 0.9, respectively. From Eq. (|4~T%|) 
we obtain W\ ~ 1 — 1.5/lng~ , W2 — 1 + 1.6/ lng -1 , 
so the width scales as W2 — w\ ~ l/lng _1 . Going back 
to unsealed velocities, the crossover takes place between 
Ci = Co — 1.5/ y/q and C2 = Co + 1.6 /y/q with a width 
c 2 - ci ~ 1/y/q. For g = 0.01, 0.001, 0.0001 one has 
c ~ 46, 218, 921 and c 2 - c\ ~ 31, 98, 310, respec- 
tively. For these high values of the velocity the distribu- 
tion function is extremely small. For instance, at c = cq, 
/(c ) ~ %exp[(q- 1 +l)]n.q + * 2 /12q-K -%\. This 
yields f(c ) - 10" 166 , 10" 2645 , lO" 36431 for q = 0.01, 
0.001, 0.0001, respectively. These values are beyond the 
accuracy of any numerical or simulation method, so the 
high energy tail in the quasi-elastic limit would look like 
a Maxwellian for the domain of velocities numerically ac- 
cessible. On the other hand, our asymptotic analysis of 
the exact solution shows that the true tail is actually 
exponential. 




0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 

w 

FIG. 3: Plot of the ratio between the velocity distribution 
function /(c) and its high energy tail Aoe _a ' c ' as a function 
of the scaled velocity w = \c\y/q/ \riq~ 1 in the quasi-elastic 
limit for q = 0.01, 0.001, and 0.0001. 



V. CONCLUSION 

The exact nonequilibrium steady state solution of 
the nonlinear Boltzmann equation for a driven one- 
dimensional inelastic Maxwell gas was obtained in 
Ref. |15| in the form of an infinite product for the Fourier 
transform (f>(k) of the distribution function /(c). The 
main goal of this paper has been to show that this rela- 
tively simple exact solution in the one-dimensional case 
also possesses the generic properties of overpopulation of 
high energy tails and exhibits a rich mathematical struc- 
ture, especially in the different limiting cases. 

We have inverted the Fourier transform to express /(c) 
in the form of an infinite series of exponentially decaying 
terms, as given by Eq. (|3.3(l with the velocity c measured 
in units of the r.m.s velocity (i.e. (c 2 ) 1 / 2 = 1). For all 
values of the coefficient of restitution < a < 1 the high 
energy tail is exponential, namely /(c) ~ Aq exp(— a\c\), 
where a = 1/y/pq = 2/ VI — a 2 and the amplitude Aq is 
given by Eq. I|3.8|l and plotted in Fig. [I] 

Special attention has been paid to two complementary 
limiting cases: the totally inelastic limit (a — > 0) and 
the quasi-elastic limit (a — > 1). In the former case some 
poles coalesce and the dominant high energy term is still 
exponential, but the subdominant term becomes an ex- 
ponential times a linear function of the velocity, where 
the numerical value of the associated amplitudes is given 
by ftlfy 

The quasi-elastic limit is much more delicate and re- 
quires some care. If we first take a — > 1 at fixed |c| 
and next \c\ — * oo (order A), the high energy tail has 
a Maxwellian form. On the other hand, if the limits are 
taken in the reverse order, i.e. first |c| — > oo at fixed a < 1 
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and then a — > 1 (order B), the asymptotic high energy 
tail is exponential. The crossover between both limiting 
behaviors is described by the coupled limit c — > oo and 
q — > with the scaling variable w = |c| ^fqj In q^ 1 — fixed 
with q = i(l — a) < 1, and occurs at w ~ 1. If id < 1 
(more specifically, 1 — w > l/lnq -1 ), the distribution 
function is essentially a Maxwellian, while the true expo- 
nential high energy tail is reached if w > 1 (more specif- 
ically, w — 1 > 1/ lnq -1 ). 

It is of interest to emphasize that the results for the 
scaling form in the quasi-elastic limit not only depend 
sensitively on the order in which both limits are taken. 
They also depend strongly on the collisional interac- 
tion, i.e. on the energy dependence of the collisional 
frequency, as well as on the mode of energy supply to 
the system. To illustrate this we have collected in Ta- 
ble Q] what is known for the different inelastic models in 
one dimension: (i) hard spheres and (ii) Maxwell mod- 
els, and for different modes of energy supply: (i) no 
energy input or free cooling, (ii) energy input or driv- 
ing through Gaussian white noise, represented by the 
forcing term —Dd 2 F(v,t)/dv 2 in the Boltzmann equa- 
tion, and (iii) energy input through a negative friction 
force oc gv/\v\, acting in the direction of the particle's 
velocity, but independent of its speed. This driving, 
referred to as gravity thermostat, can be represented 
as the forcing term g(d/dv)[(v/\v\)F(v, t)] in the Boltz- 
mann equation. The results corresponding to order A 
with the gravity thermostat have been obtained by the 
same method as followed in Ref. 8]. It is worthwhile 
noting that in the quasi-elastic limit a bimodal distri- 
bution, i [5(c + 1) + S(c — 1)], is observed in inelastic 
hard sphere systems both for free cooling and for driv- 
ing through the gravity thermostat, whereas in inelas- 
tic Maxwell models this bimodal distribution is only ob- 
served for the gravity thermostat. 

It is important to note that in the normalization where 
velocities are measured in units of the r.m.s. velocity, the 
high energy tail in the driven inelastic Maxwell model is 
only observable for very large velocities, as illustrated in 
Figure |21 for strong (a — > 0) and intermediate (a = i) 
inelasticity. In the quasi-elastic limit, where (a — > 1), the 
tail is even pushed further out towards infinity, as ana- 
lyzed at the end of Section IIVI This also explains how 
to reconcile the paradoxical results of exponential large- 
c behavior with the very accurate representation (|4.4() 
of the distribution function in the thermal range, in the 
form of a Maxwellian, multiplied by a polynomial expan- 
sion in Hermite or Sonine polynomials with coefficients 
related to the cumulants. The validity of these polyno- 
mial expansions, over a large range of inelasticities with 
(0 < a < 1) had been observed before in Q for inelastic 
hard spheres and in 01 f° r inelastic Maxwell models. On 
the other hand the high energy tail is cx e~ a l c l, and not 
cx c N e~~ c / 2 , where N is some large number, and yields 
diverging moments M 2 „ = (c 2n ) and cumulants Ca« in 
the limit n — > oo, as shown in Section ITT1 

The exact solutions of the nonlinear Boltzmann equa- 



tion for the freely evolving |13| and the driven [15J in- 
elastic Maxwell model (extended in this paper), as well 
as the rigorous proof of [24j f° r the long time approach 
of the distribution function to a scaling form, validate 
the self consistent method, developed in |5( for analyti- 
cal studies of possible over- or underpopulations of the 
high energy tail of velocity distributions, not only for in- 
elastic Maxwell models, but - more importantly - also 
for inelastic hard sphere, where exact solutions are not 
known. This possibility of assessing the validity of gen- 
eral kinetic theory methods by means of exact solutions 
of the nonlinear Boltzmann equation is one of the main 
reasons why the study of inelastic Maxwell models is of 
interest. 
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APPENDIX A: LARGE-fc EXPANSION 

The asymptotic behavior of ip for large k can be ob- 
tained by inserting the ansatz, ip = — X\k\ + ln(Ak 2 ) + 
Y^=i a nk~ 2n , with unknown coefficients {A,a n } into 
Q2.2|) . and equating the coefficients of equal powers of 
In k and k n with the result, 

m = -x\ k \ + H k 2 / Pq ) - £ { ~l ) \^ + pq J2_ v 

n—l 

( A1 ) 

where A is as yet undetermined. The series converges 
for k 2 > q/p. However for the tp(k) above to qualify as 
a solution of (|2.2I) the radius of convergence is further 
restricted to (qk) 2 > q/p or pqk 2 > 1. The constant 
A must be chosen such that ip(k) satisfies the boundary 
condition ip ~ ™i^ 2 a t small k. This can be done by 
matching QAlfl with 12.4fl . The latter satisfies already 
the small-fc boundary condition. Matching at pqk 2 = 1 
yields then, 

^= -21n(p,)+E~ 1 ^ 1 
x (j-pk^iK + p -2 n+ 1 q -2 n _ 1 'j ■ (A2) 

Both terms can be combined into a single n-sum with 
n = ±1, ±2, • • • . The above result is not only convenient 
for numerical evaluation, as shown in Fig. 0] but also 
for analytic evaluation in two limiting cases. We first 
consider the totally inelastic limit (a — > or p = q = ^). 
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TABLE I: Asymptotic behavior of the distribution function /(c) for one- dimensional systems in the quasi-elastic limit. In 
general, the result depends on the order of limits. Order A corresponds to take first a — > 1~ and then cj — > oo, whereas order 
B refers to the reverse order, i.e. first jcj — *■ oo and then a — *■ 1~. The first/second footnote in the second column gives the 
reference where the result for order A/B was obtained. 



State 



System 



Order A 



Order B 



Free cooling 
White noise 
Gravity thermostat 



a Ref. 
<>Ref. 
c Refs. |] 
d Refs 
This work 
■friefs. IIE 
9Ref. [ia 

h Ref. Ea 




Hard spheres" 6 
Maxwell model 00 
Hard spheres ad 
Maxwell model 6 ^ 
Hard spheres' 39 
Maxwell model 6 ' 1 



i[5(c-l)+(5(c+l)] 



e 
e~ 



I[i(c-l)+tf(c+l)] 
i[5( C -l) + 5( C +l)] 



c 

-a|c| 3 / 2 
,,,-'> I <■ I 



-a|c| 




a 



FIG. 4: Coefficient A as a function of the coefficient of resti- 
tution. The arrow indicates the value A = 7r/21n2 at a = 0. 
The dotted line represents the asymptotic form Eq. 1A7I for 
small q. 



There the expansion IjAlfl can be cast into a simpler form, 



-Mk\ + ln(4fc 2 ) - |S^i ( "„ r i-2-('" + i) 
-A|fe| + ln(4fc 2 ) + 1E™ =0 2^ ln(l + ^) 
-A|^| + iE™ = o^ln(l + 2 2m fc 2 ). (A3) 



Matching this expression in k 2 = 1 /pq = 4 with the exact 
solution in (£3, V'C 2 ) = ~ Em=o 2™ ln(l + 2- 2m ), yields 



the nice result 
A 



1 

41n2 



m— — oo 
oo 



2 m ln(l + 2- 2m ) 
dxx~ 3/2 ln(l + or) = 



21n2 



(A4) 



One can verify using the Euler-MacLaurin summation 
formula (see Eq. (23.1.30) of j2^) that all correction 
terms to the integral are vanishing, and the integral is 
listed in Eq. (4.293.3) of j^. 

In the quasi-elastic limit (a — > 1 or q — > 0) the sum 
originating from the second term inside (• • • ) in ljA2fl is of 
0(q 2 ), and will be neglected. To evaluate the first term 
T for small q we expand it as follows, 



E 

fc=i 
00 

E 



(-1) 



k+l 



1 



1 — p 



:2 k 



■2 k 



(-1) 



fe+i r 



1 

~2~kq 



1 4 

2/c 



1 — j) 



2fr 



2(2fc-l) 



■2k 



4k 



{2kqf 



0(q) 



2q 



Ii 2 (-l)-±Lii(-l) + ±Li 2 (-l) 



0(q), 
(A5) 



where the poly logarithmic functions are defined as, 



Li fc (z) = 5>"A 



with Li 2 (-1) = -j^tt 2 and Li x (-1) = - la 2 H3 
final result for A at small q is then, 



(A6) 
The 



A = 



1 



24 



-2qlnq+± 9 (ln2+±- f|) 



+9 2 ln g + C( g 2 



(A7) 
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APPENDIX B: ASYMPTOTICS IN 
QUASI-ELASTIC LIMIT 

To calculate in (|4.7fl for small q we expand the log- 
arithm, and perform the (m, ^-summation. The result 
is, 



The small- x expansion of l/(e x — 1) has been obtained 
from (|B4[1 with x — > — x. 

By combining the relations (|B6|) and (|B7|I with the 
small- a: expansion of ln(l — e~ x ), we obtain from l|B5(l . 



r(3) _ v^oo 

D » ~ l^k=i k(i 



v 



2k,„2ksn+l 



_+il n( l_ e -x ) _I Li2(1 _ e -x ) 



E 



fe=l fc(l- p 2fc) 



, _l ("+l)g 2 



S(a:)+<JS(n,a;)+o(g), 



(Bl) 



E 



£?2r. 



^ 2m(2m- 1) 



(1 - hx5 m l) 



where all contributions oc q° and oc q have been included. 
The dominant term is, 



= — + ±\nx-l + R -£ + O(x 2 ), (B8) 



S(x) 



OC 

E 

fc 



' k(l - e~ kx ) 



(x=-21np). (B2) 



where R = £~ =1 B 2m /[2m(2m - 1)]. As \B 2k \ ~ 
2(2fc]!/(27r) 2fe for fc -> oo (seeEqs. (23.2.16) and (23.2.18) 
of 29] ), the series is a divergent asymptotic series with 
alternating signs. One obtains the greatest accuracy, 
denoted by i? mo) , if one breaks off the series just be- 
fore the smallest term in the series, which is defined to 
be the (mo + l)-th term. Then the maximum error is 
SS{x, n) = \ + \q(n + §) ~ \ + \x{n + §). (B3) \B 2mo+2 / (2m + l)(2m + 2))| (see Chap. 0.33 of 

In the present case one can simply verify that mo = 3, 
and the best possible estimate for the remainder in the 

kx \ 

) limit where q — > is given by, 



In the remaining contributions to (|B1|I only the term k 
1 needs to be taken into account, and yields 



To study the small- a; behavior of (|B2Jl we construct an 
asymptotic series for S(x), by expanding 1/(1 - 
in powers of x. This can be done most conveniently by 
using the small- a: expansion of xcothx, or equivalcntly 
the generating function for the Bernoulli numbers, B 2k = 
(~l) k+1 \B 2k \ (see Eqs (23.1.1-2) of [U), which we write 
as, 



1 



1 



1 

x 



E 



B 



2m x 2m-l 



Substitution of (|B4(I with x 
1 



S(x) 



= x ( 2m )- 
kx into (|B2|) yields, 
1 



(B4) 



Li 2 (e- X ) - -ln(l - e- x ) 



00 T 2m-1 

(2m)! 

m— 1 v 7 



-S2mLi2-2m(e 



(B5) 



We have used the definition i|A6(l of the polylogarithmic 
functions, which are all singular in x = 0. To deter- 
mine the behavior of the dilogarithm ; L i 2 (e~ x ) we use 
the functional relation (see Eq. (5) of |31|). 



0.0813492 ± 0.0005952. (B9) 



The inaccuracy in S(x), caused by the inaccuracy in the 
asymptotic series i?0i can be substantially reduced - if 
so desired - by restricting the m-sum in (|B5ll to mo = 3 
terms, and calculating the difference, 



A(a:) = £ 



k=l 



l-e- 



1 1 

kx ~ 2 



E ^ kx ) 

m=l y ' 



(B10) 



in the small- a; limit as an integral. The result for in- 
stance, to seven decimal points is A = —0.0002877. 
Hence 



Li 2 ( e - X ) = \ - ln(e- x ) ln(l - e~ x ) - Li 2 (l - e~ x ) 
= ^. + x (lnx- ±x) - (x - \x 2 ) + 0(x 3 ). (B6) 



Here the small- a; expansion of the sum in i|B5l) can be 
obtained from the relation, 



Ro = i^ 3) + A ~ 0.0810615. 



(Bll) 



Combination of the results I|B1J1 , (|B3|) and l|B8(l gives the 

(3) 

dominant small- a; behavior of _B„ in the form, 



Li_ n (e-) = Er=ifc T ' 



_ (_d_\ n f 1 _ 1 , v^oo B 2m „2m-ll 
— V dx) \ x 2 ^ Z^m=l (2m)! X J 

= [l - \x8na + o(x)] . 



(B7) 



B { V = ^ + \\nx-l + R + \x(n+^). (B12) 

Final elimination of x = 2q(l + ^q + ^q 2 + ■ ■ ■ ) in favor 
of q gives (|4.9|) in the main text. 
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